Nanoscale structures formed in silicon cleavage studied with large-scale electronic 
structure calculations; surface reconstruction, step and bending. 
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The 10-nm-scale structure in silicon cleavage is studied by the quantum mechanical calculations 
for large-scale electronic structure. The cleavage process on the order of 10 ps shows surface re- 
construction and step formation. These processes are studied by analyzing electronic freedom and 
compared with STM experiments. The discussion presents the stability mechanism of the experimen- 
tally observed mode, the (lll)-(2 x 1) mode, beyond the traditional approach with surface energy. 
Moreover, in several results, the cleavage path is bent into the experimentally observed planes, owing 
to the relative stability among different cleavage modes. Finally, several common aspects between 
cleavage and other phenomena are discussed from the viewpoints of the nonequilibrium process and 
the 10-nm-scale structure. 

PACS numbers: 71.15.Pd, 68.35.-p, 62.20.Mk 



I. INTRODUCTION 

Cleavage is a nonequilibrium process and its dynami- 
cal mechanism is essential. In particular, the cleavage of 
silicon single crystals is of great interest from the multi- 
scale viewpoint between macroscale and atomicscale pic- 
tures. In the macroscale picture, silicon shows perfect 
brittleness and brittle fracture is usually described by 
continuum mechanicsi^ In the atomicscale picture, on 
the other hand, a cleaved surface contains areas with 
well-defined reconstructions. Currently, these atomic- 
scale structures are observed by scanning tunneling mi- 
croscopy (STM)4 and other experiments. The multiscale 
feature of the phenomenon, as discussed in this paper, 
appears in nanoscale (or 10 nm scale) processes near the 
crack tip, though such processes cannot be seen by direct 
(in situ) experimental observation. 

A fundamental question is what Miller index and sur- 
face reconstruction appear at the cleavage surface. A 
traditional prediction is that the cleavage plane should 
be that with the smallest surface energy, or the small- 
est energy loss with forming surface. This prediction, 
however, is actually not satisfactory owing to the follow- 
ing experimental facts; (i) The easiest cleaved surface 
of Si is a metastable (lll)-(2 x 1) structuro 4 i 5 i 6 i 7 i 8 i 9 i 10 
and will change, irreversibly, to the ground-state (7 x 7) 
structureAii (ii) The (110) cleavage plane is also experi- 
mentally observed but less favorable (see an experimental 
(STM) study, — and a recent theoretical studj*i2>). (iii) 
The cleaved Ge(lll) surface also shows the same (2x1) 
structure, while the ground state surface structure is the 
c(2 x 8) structurei 4 * 1 ^ 1 ^ These facts imply the importance 
of direct cleavage simulations with electronic structure. 
Although such simulations have been carried out thus 
£ ar( 8 i 9 i i3 |- ne investigation is still limited, owing to the 
system size of 10 2 atoms. 

In this paper, the cleavage of silicon is studied with 
quantum mechanical calculations for large-scale elec- 



tronic structuresi£iiLi2ii2i22i2i and we use a transferable 
Hamiltonian 2 ^ in the Slater-Koster (tight-binding) form. 
The methodology is reviewed briefly in Appendix^ The 
method realizes the cleavage process with more than 10 4 
atoms or a sample length of 10 nm. This paper is orga- 
nized as follows; Section [n] describes the important as- 
pects discussed in this paper. The easiest cleavage mode 
on the (lll)-(2x 1) plane is discussed in Section ITTT1 and 
Section llVl The latter section focuses on step formations. 
The simulation results for bending in the cleavage path 
are presented in Section^ Finally, in Sect ion l*VTl several 
common aspects of nanoscale structures are discussed for 
cleavage and other phenomena. 
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FIG. 1: (a)-(c): Toy model of cleavage, in which surface 
reconstruction is illustrated as dimerization. (d): Possible 
cleavage paths on Si(lll) plane. 
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FIG. 2: Cleavage process of silicon with (lll)-(2 x 1) cleaved surface. (a)-(e):Successive snapshots with a time interval of 
approximately 2 ps. (f): A part of the lower cleavage surface shown in the snapshot (e). A step is found and is classified into 
the l [211]-type' or 'via-(lll)-plane type' that can be decomposed into the successive bending of cleavage planes as (111) — > 
(111) -(111). 



II. ASPECTS OF CLEAVAGE PROCESS IN 
REAL CRYSTAL 

The first aspect of cleavage is the typical time scale as a 
nonequilibrium process. The time scale is determined by 
the cleavage propagation velocity. In the continuum me- 
chanics and many experiments, the propagation velocity 
is given on the order of, but less than, the sound velocity 
or the Rayleigh wave velocity (wr = 4.5km/s = 4.5nm/ps 
for Si) ^iiSi Since the atomistic process occurs within the 
time scale, the reconstruction in cleavage, unlike that 
in annealing, should occur locally. In a typical pro- 
cess, a surface-bound state is formed between electrons 
in nearest-neighbor dangling bond sites, as illustrated in 
Figs. IHa)-(c). In other words, the elementary process 
should contain two successive bond breakings, not one. 
This nearest-neighbor reconstruction mechanism will di- 
rectly give the experimentally observed (2 x 1) recon- 
struction (See Section ITTTf> . 

The second aspect comes from the crystal structure. 
Figure IHd) shows possible cleavage paths on the Si(lll) 
plane, with or without step formation. Unlike the toy 
model of Figs. [IJa)-(c), the system does not have the 
mirror symmetry with the cleavage plane and the upper 
and lower cleaved surfaces are inequivalent in symmetry. 
In particular, the inequivalence between the upper and 
lower stepped paths in Fig. ^d) is distinctive in experi- 
ments. 

Since experiments reported only the (111) or (110) 
cleavage plane, theory should explain why other sur- 
faces do not appear. An interesting fact is that the cal- 
culated surface energy of the reconstructed (001) sur- 
face is smaller than that of the (lll)-(2 x 1) surface 
(T^) = l-41J/m 2 < 7^ = l^J/m 2 )^. This 
means that the absence of the (001) cleavage surface is 



not simply predicted by surface energy^ The possibil- 
ity of the (001) cleavage mode was investigated in our 
previous studyr& in which the (001) cleavage planes are 
observed for small sample sizes less than 10 nm. For 
larger sizes, the flat (001) cleavage surface becomes fairly 
unstable and many steps form. The instability of the 
(001) cleavage mode is helpful in understanding the sta- 
bility of the (111) cleavage mode, as discussed in Section 
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III. (lll)-(2 x 1) CLEAVAGE MODE 

In experiments, the cleavage on the Si(lll)-(2 x 1) 
plane is the easiest cleavage mode. As the atomic struc- 
ture of the cleaved surface, Pandey's 7r-bonded structure 5 
is now widely acceptedA&L&Siiii An actual surface for- 
mation process was shown^ when a parallel separation 
is introduced in a slab with the minimal periodic simu- 
lation cell for the (2 x 1) structure (See Figs. 1 and 2 of 
Refi&). 

The present cleavage simulations are performed with 
the [211] propagation directions (Fig. ^d) for geome- 
try), which are consistent with those in typical experi- 
ments (See Fig. 1 of Ref. 26 , for example). An external 
load is imposed and its physical origin is the concen- 
trated elastic field in the macroscale experimental sam- 
ple. See Appendix iBl for details. A smaller sample with 
416 atoms was also simulated as shown in Appendix IU1 
so as to confirm the quantitative agreement with pre- 
vious experimental or theoretical studies. In all cases, 
the periodic boundary condition is imposed in the [011] 
direction, which is orthogonal to the cleavage propaga- 
tion direction. The periodic length contains eight atomic 
layers or is four times larger than that of the minimum 
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FIG. 3: (a)-(d):Initiation of cleavage while forming the (lll)-(2 x 1) structure, (e)-(f): Schematic figures of the transition 
from (e) the buckled (2 x 1) structure to (f) the (tilted) 7r-bonded (2 x 1) structure. The former structure appears on the 
lower cleaved surface in (b) and the latter appears in (d). The oval in (e) indicates the presence of a lone pair state on the D 
atom site. The oval in (f) indicates the 7r-bonding zigzag chain in the direction perpendicular to the page, (g)-(i): Quantum 
mechanical analysis of a process in which the bonding wavefunction (</>;) between A and B sites changes into that between A 
and C sites; (g) The one-electron energy Si = (cf>i\H\cf)i) and the weight of s orbitals fg for the wavefunction <j>i. (h) The 
atomic distance between A and B sites (cIab) and that between A and C sites (cLac)- (i) The spatial weight distribution on 
the A, B and C atom sites (riA, ub, nc) for the wavefunction <f>i. The rest of the weight (1 — n\ — hb — nc) is also plotted. 



crystalline periodicity. In some simulations, the atomic 
motion is under a constraint by the minimum crystalline 
periodicity in the [Oil] direction. We call this constraint, 
the '2D-like constraint'. For a systematic investigation, 
simulations were carried out with and without the '2D- 
like' constraint. 

Figurc|21shows an example of the simulation within the 
'2D-likc' constraint, in which a step appears. The sam- 
ple contains 11,096 atoms. Hereafter, the bonds (rods) 
are drawn just as a guide for eye. From Fig. [21 the 
cleavage propagation velocity is estimated to be v plop ~ 
2nm/ps=2km/s, as expected in Section^] The cleaved 
surface in Fig.^f) contains the 7r-bonded (lll)-(2 x 1) 
structure, of which unit structure is a set of seven- and 
five-membered rings. 

A. Elementary reconstruction process 

Figures|3Ja)-(d) show an example of (2 x 1) reconstruc- 
tion. As in the previous simulation with the parallel 
separation^ the reconstruction process has two stages; 
The buckled (2 x 1) structure of Fig. |3{e) appears, as 
shown in Fig. E^b), when the two dangling-bond elec- 
trons at the C and D atom sites form an atomic (lone- 
pair) state at the D atom site^ After that, the 7r-bonded 
(2 x 1) structure of Fig- Etf ) appears, since (i) the atomic 
state is transformed into the 7r-bonding state between the 



B and D sites and (ii) the bonding state between the A 
and B sites is transformed into one between the A and C 
sites, so as to reproduce the tetrahedral coordination in 
the C site. In the three-dimensional view, a 7r-bonding 
zigzag chain appears in the direction perpendicular to the 
page. Hereafter, the transition from the buckled struc- 
ture into the 7r-bonded one is called the 'BP transition'. 
Note that the resultant 7r-bonded structures are always 
tilted, since the D or B atom site is shifted to the vacuum 
side. The two tilted structures arc incquivalent in sym- 
metry and have a small difference in surface energy^fl^ 
In this paper, however, we do not focus on the difference 
between the two tilted structures. 

The above explanation of (ii) is confirmed by the fol- 
lowing quantum mechanical analysis; Figure EUg) shows 
the one-electron energy ei = (4>i\H\<fri) and the weight 

of s orbitals For example, fs^ — 1/4 in an ideal 

sp 3 hybridized state. Figure BIh) shows the distances 
between sites A, B and C. Figure [3{i) shows the spa- 
tial weight distribution (|0i(t")| 2 ) on atom sites. Note 
that the intermediate wavefunction, indicated as the ar- 
row in Fig.[3^h) or (i), shows no characteristic feature in 

(i) 

£i and fs , such as maxima, minima or plateaus. This 
is quite different from the reconstruction process on the 
(001) surface^ in which the intermediate wavefunction 
shows a plateau with a large weight on the s orbitals 
« 0.8). 
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FIG. 4: (a) Si(lll) surface with coexistence of the buck- 
led and 7r-bonded (2 x 1) structures. The induced strain is 
shown as an arrow at the boundary between the buckled and 
7r-bonded structures. The dashed line indicates the bond site 
that will appear after the next induced BP transition, (b)- 
(c): Two snapshots of cleavage process. The solid arrow cor- 
responds to the direction of the possible surface strain force 
shown in (a). The dashed arrow denotes the cleavage propa- 
gation direction. 



B. Role of anisotropic surface strains 

We studied the role of anisotropic surface strains in the 
reconstruction process. As a demonstration, Fig. 0{a) 
shows the (111) surface with the coexistence of the buck- 
led and 7r-bonded (2x1) structures. An anisotropic strain 
force is induced, which is depicted by an arrow. The 
strain force enhances a further BP transition of the next 
left unit. Owing to the crystalline symmetry, the strain 
force on the upper and lower cleaved surfaces shows the 
opposite direction to each other, as shown by solid ar- 
rows in Fig.QJb). This explains that the 7r-bonded struc- 
ture appears on the lower cleavage surface in Fig. EJd), 
but not on the upper one. Since the experiments on 
cleaved samples support the 7r-bonded structure^iiSi the 
7r-bonded structure should appear commonly on the up- 
per and lower surfaces. Such a result is obtained in 
Fig. EKc), m which an additional constant force is im- 
posed in the left direction on several leftmost atoms of 
the upper cleaved surface. The additional force is created 
by the boundary condition of surface strain, since the 
simulation sample should be embedded in a macroscale 
sample. In the embedded situation, the boundary region 
of the simulation sample should be connected with suc- 
cessive 7r-bonded (2x1) structures and should be under 
the resultant strain force. In Fig.^Jc), the resultant sam- 
ple contains the 7r-bonded (2x1) structures on the upper 



and lower surfaces. 

We also found that the sample geometry can have an 
effect on the possible BP transition, through strain. For 
example, the BP transition seems to occur more easily 
in a thicker sample, with a larger sample length in the 
[111] direction, since the BP transition should accompany 
the strain in deeper (subsurface) layers. This observation 
was unchanged, when additional simulations were carried 
out by more complicated Hamiltonians with a better re- 
production of the surface energy. A more systematic in- 
vestigation should be carried out, in the future, on the 
quantitative condition of the BP transition. 



IV. STEP FORMATION IN CLEAVAGE 

Several (111) cleavage simulations contain step forma- 
tion. There are two inequivalent paths shown as the up- 
per and lower paths in Fig.^d). The 'upper' type of path 
is usually called the '[211] type', since the step is descend- 
ing in the [211] direction. In the same manner, the 'lower' 
type is called the '[211] type'. The two types of step are 
experimentally reported. In early papers^*^*^ it was re- 
ported that the lower path in Fig.QId), the '[211] type' 
is predominant. In later experimental papers^, however, 
the upper path, i.e. the '[211]' type step, was reported. 
In our simulations, both types of step are observed. 



A. 'Upper'-type step 

First, we discuss the step in Fig. [21 which is classified 
into the 'upper' path in Fig.^d) or the [211] type. This 
type of step was observed with and without the '2D-likc' 
constraint. The formation process is shown in Fig. [3] 
Hereafter, the members of several rings are plotted, such 
as '6', '5' or '7', so as to clarify the reconstruction. This 
step can be decomposed into the successive bendings in 
the cleavage path as (111) — > (111) — > (111) planes (See 
Fig.[IJd) and Fig.(2J|. That is, the step is formed through 
the (111) plane and we call the step, the 'via-(lll)-plane' 
type. In the crystalline geometry, the (111) and (111) 
planes are equivalent. If we ignore the reconstruction (or 
quantum mechanical) freedom, the step formation can 
be understood, since bending between equivalent planes 
can occur by the local fluctuation of the (concentrated) 
elastic field at the crack tip, particularly in its angular 
dependence. 

As a quantum-mechanical analysis, snapshots of the 
step formation are shown in Fig. [5] in which an atom 
with an excess or deficit electron population is marked 
'+' or '-', respectively. The initial bonding states (dashed 
line) changed into the surface-bound state mainly local- 
ized on a '+' atom site. An atom '+' is always placed at 
the vacuum side of a buckled structure, which can be un- 
derstood by the general quantum mechanical tendency. 30 
The resultant surface shows a balanced structure, owing 
to the alternate alignment of the '+' and '-' sites. The 
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present step structure is concluded to be a possible one, 
although it has not yet been confirmed experimentally. 

Now we comment on a pioneering theoretical paper, 
by Chadi and Chelikowsky,— in which the 'upper' and 
'lower' types of step are compared with assumed atomic 
structures. They concluded that the upper path is unreal- 
istic, which seems to be in contrast to the present result. 
The contrast appears, because the present step structure 
is different from that assumed in Refi^i In this reference, 
the dangling bonds at the step edge are assumed to be 
rebonded with each other in the perpendicular direction 
of the page of Fig. IHd), but such a rebonding process 
results in only a tiny energy gain. In the present step 
structure of Fig. [SJ on the other hand, the dangling bond 
electrons at the step edge are rebonded within the plane 
of Fig. Did), whose energy gain mechanism is explained 
above. 





FIG. 5: Step formation process on (111) cleavage plane. The 
step is classified into the '[211]-' or 'via-(lll)-plane'-type step 
(See Fie.fflf)). Two snapshots with a time interval of approx- 
imately 0.6 ps are shown. The dashed lines indicate the initial 
(crystalline) bonds. The '+' or '-' symbols on atoms indicate 
the excess and deficit electron populations, respectively. 



B. 'Lower '-type step 

Several simulations result in the appearance of the 
other type of step, the lower-type step in Fig. ^d) or 
the '[211]' type, when they were carried out without the 
'2D-like' constraint or minimal periodicity in the [Oil] di- 
rection. An example is shown in Fig. EI a). To observe the 
freedom in the [Oil] direction, we classify the atoms, by 
their initial positions, into subsystems or 'slices'. Each 
'slice' contains atoms within the minimum unit length 
(two atom layers) in the [Oil] direction. Four slices 
are dchncd in Fig. EJa) and two of them are shown in 
Figs.OJb) and (c). 

In Fig.UJa), the cleaved surface is defective before the 
step formation. In the defective area of Figs.|£Jb) and (c), 
six-membered rings appear in the surface layer, as well 
as five- or seven-membered rings. A defective area also 
appears in flat (nonstepped) areas, as in Fig. The ap- 
pearance of the defective six-membered rings means that 
the reconstruction of the dangling-bond electrons does 
not occur within the slice. We should recall that an ideal 
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FIG. 6: (a) (111) cleaved surface with '[211]'- or 'via-(lOO)- 
plane'-type step. The arrow indicates the cleavage propaga- 
tion direction. (b)(c): Two 'slices' of (a) (See text for details). 



(Ill) surface has symmetry with a ±27r/3 rotation and 
the (2 x 1) reconstruction mechanism shown in FigUfe) 
or (f) is possible in the [211] , [121] and [112] directions. 

This type of step formation is decomposed into the suc- 
cessive bendings of cleavage as (111) — > (100) — > (111) 
planes. In other words, the step is formed through the 
(100) plane and we call the step the 'via-(100)-plane' type 
(See Fig.^d) for geometry). At the step edge, the dan- 
gling bond sites are equivalent to those on a (100) or (001) 
plane. The reconstruction among them is possible with 
an energy gain of dimerization^l The actual process of 
successive bond breakings and reconstructions on (100) 
surface can be obtained as an unstable cleavage mode^ 
With the 2D-like constraint, however, the dimerization is 
prohibited. This is why this type of step does not appear 
within the 2D-like constraint. 

The geometries in Fig. can be compared with 
those proposed in STM experiments, 28 The geometry 
of Fig. Etb) is identical to that of Fig. 5(b) in Ref. 2 *, 
because of the same alignment of the ring structures; 
5 — > 7 — * 6 — * 5 before the step and 6 —> 5 — > 7 after 
the step. In the same manner, the geometry of Fig. (Bfc) 
is identical to that of Fig. 5(a) in Refi2£ only before the 
step (5 — * 7 — > 5). After the step, however, the present 
geometry (6 — »• 5 — > 7) is different, by one six-membered 
ring, from that shown in Ref^& (6 — > 6 — > 5 — > 7). 

The present investigation does not show that the 
upper- or lower- type of step in Fig.^d) is more favorable 
than the other, because the present simulations were car- 
ried out under the periodic boundary conditions in the 
[011] direction. We speculate that the appearance of the 
upper path, the l [211]-type' step, should be enhanced, 
because the step can be formed within the minimal crys- 
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FIG. 7: Defective area on cleaved Si(lll) surface. Three 
slices are picked out. The six-membered rings marked (I) and 
(II) are geometrically equivalent among slices (a)-(c). 

talline periodicity in the [Oil] direction. We have been 
informed that the two types of step are observed in STM 
images^ and we consider that a more systematic investi- 
gation is desirable both theoretically and experimentally. 

C. Step formation and stability of cleavage mode 

The present study of step formation can lead to the 
stability mechanisms of the (lll)-(2 x 1) cleavage mode, 
since a stable cleavage mode should be robust against pos- 
sible disorders and fluctuations. We have discussed that 
a step formation is decomposed into successive bendings 
of cleavage plane, such as (111) — > (100) — > (111) planes. 
The former and latter bending processes correspond to 
the deviation and (quick) recovery of the (111) cleavage 
mode, respectively. Since the 7r-bonded (2x1) structure 
has a 7r-bonded zig-zag chain perpendicular to the cleav- 
age propagation direction, it can accompany the order- 
ing in the surface structure. Actually, in Fig. [5] the step 
formation changes the cleaved surface from a disordered 
(defective) structure into an ordered (2 x 1) structure. 

Another stability mechanism can be found in the pos- 
sibility of multiple propagation directions. When the 
(2 x 1) structures on the (111) and (001) surfaces are 
compared, a crucial difference comes from their symme- 
try. In the limited growth of the (001) cleavage mode^ 
the cleavage propagates easily in the dimerization direc- 
tion. The axis of the dimerization direction is unique on 
the (001) surface and a (single) step formation changes 
the axis perpendicularly ((2x1)— > (1 x 2)). As a conse- 



quence, a cleavage growth with multiple propagation di- 
rections requires step formation and the resultant (001) 
cleaved surfaces tend to be rough with many steps^ The 
steps on the (111) surface, on the other hand, can pre- 
serve the cleavage propagation direction or the direction 
of the 7r-bonded chain ((2 x 1) — > (2 x 1)), as shown in 
this section. Moreover, the ideal (111) surface, unlike the 
(001) surface, has three equivalent axes of stable cleavage 
propagation; the [211], [121] and [112] axes. In conclu- 
sion, the stability of the (111) cleavage mode is supported 
by the geometrical features; (i) A step can be formed 
without changing the cleavage propagation direction, (ii) 
The cleavage propagation direction can be changed with- 
out step formation. The above feature (ii) should result 
in multiple domains without step formation. Actually, 
an experimental STM image, Fig. 2 on Refi2&, shows a 
cleaved (111) surface with multiple domains. We spec- 
ulate that such a multiple domain structure is formed, 
because the cleavage propagation directions are locally 
different. 

We should emphasize that the present stability mecha- 
nisms, the robustness against step formation and the pos- 
sible growth with multiple propagation directions, cannot 
be explained by traditional discussion on surface energy. 



V. BENDING IN CLEAVAGE PATH 

Finally, a direct investigation was carried out for the 
preference of the (111) or (110) surface in the cleavage 
process. We expect that, even if the simulation initiates 
the cleavage on other planes, the cleavage path will be 
bent into a more favorable plane. The above was con- 
firmed by simulations. Hereafter, symmetrically equiva- 
lent planes are called by their 'type'. For example, (111) 
and (111) planes are (lll)-type planes. 

An example is shown in Figs. IHIa)-(e). The sample 
contains 10368 atoms and the sample size is lOnmx 14nm 
on the (110) plane. The periodic boundary is imposed in 
the [110] direction by eight atomic layers. We prepare an 
initial 'seed' of cleavage on the (001) plane. An external 
force acts on the atoms near the sample surfaces except 
the left one. See Appendix iBl for details. 

Figure OJf) shows the broken-bond sites on the crys- 
talline geometry and the corresponding indices are drawn 
in Fig.|Sfg). In this case, the bendings are observed from 
the initial (001) plane to (lll)-type and/or (llO)-type 
planes. The actually observed path shows the bending 
of (001) -> (111) -> (110) planes. Here, the (110)- 
type planes appear only near the sample surface, which 
implies that the appearance of (llO)-type planes is en- 
hanced by the present sample geometry. Different paths 
were observed under different simulation conditions, such 
as tuning the magnitude of the external load, preparing 
different sample geometry, changing the region of impos- 
ing external force and setting a different cleavage seed. 
For example, we observed a bending between two (111)- 
type planes, such as (111) — > (111) planes. In conclu- 
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FIG. 8: Cleavage process of silicon with bending of cleavage planes, (a)-(e): Snapshots with a time interval of approximately 1 
ps. Atoms are drawn as balls from the viewpoint of the [110] direction, (f) : The broken bond sites in snapshot (e) are plotted 
as bold lines in the ideal (crystalline) geometry, (g) Indices of geometry in (f ) . 



sion, a (lll)-type or (llO)-type plane appears, while no 
(OOl)-type plane appears except the initially prepared 
one. This conclusion is consistent with the experimental 
preference of the (111)- and (llO)-type cleaved surfaces. 

Moreover, the bending path changes, sometimes, from 
a disordered (defective) area into an ordered one with a 
well-defined reconstruction; In Fig. 02a), the (110) plane 
after bending contains the buckled structures, known 
for typical reconstruction^* In the process shown in 
Figs.|5fb)-(c), the resultant (111) cleavage plane contains 
the 7r-bonded (lll)-(2 x 1) structure. 



VI. GENERAL DISCUSSION 

The present investigation of cleavage can be discussed 
from the general viewpoints of (I) time scale and (II) 
length scale. From the viewpoint of time scale (I), the 
metastable Si(lll)-(2 x 1) structure appears, owing to 
the local reconstruction mechanism within a limited time 
scale. Since the mechanism is general in the quantum me- 
chanical picture, we expect the (lll)-(2 x 1) surface also 
in other noncquilibrium processes. Actually, we were in- 
formed that the corresponding STM image is found in 
a Si(lll)-\/3 x y/3-Ga. surface, when the Ga atoms are 
locally removed at room temperature.™ At higher tem- 
peratures, the (7 x 7) structure appears^ The above sit- 
uation is another 'fast' process almost free from thermal 
equilibration. A review article 4 lists other experiments 
with the appearance of the Si(lll)-(2 x 1) structure. 

From the viewpoint of length scale (II), structure of 
material is generally determined by volume and sur- 
face terms of energy. Dimensional analysis shows a 
crossover of mechanical property between nano- and 
macroscale systems; The two energy terms are competi- 
tive in nanoscale system, while the volume term is dom- 
inant in macroscale system. In the case of silicon nanos- 
tructure, the above energy competition originates from 
that between the sp 3 (bulk) term and the non-sp 3 (sur- 



face) term. In the present investigations, the crossover 
at the 10 nm scale is found in step formation or bend- 
ing in cleavage path. We speculate that the crossover 
at the 10 nm scale is universal in silicon, because the 
surface energy is always on the same order among differ- 
ent surface indices and reconstructions^* The crossover 
may be seen also in the size dependence of the shape 
in self-organized Si islands^* Islands with a size of 10 
nm or less have a semispherical shape and those with a 
larger size have a pyramidal shape with facets in well- 
defined indices. These crossovers at the 10 nm scale can 
be understood, because a well-defined reconstructed sur- 
face appears, only when the system size is much larger 
than the unit length of the reconstructed surface (~1 
nm). 

The present method of large-scale electronic structure 
calculation has wide applications and is not specific to 
cleavage. The structural property at the 10 nm scale is 
an urgent problem of the present technology and can be 
a future study of the present method. 



APPENDIX A: METHODS FOR LARGE-SCALE 
ELECTRONIC STRUCTURE CALCULATION 

In recent years, we have developed theories and 
program codes for large-scale electronic structure 
calculations>i§iiLi2ii2i22i2i Among them, a variational 
(VR) procedure for generalized Wannier state d 16 ! 21 is 
used in the present cleavage simulations. The Wannier 
state is well-defined localized 'chemical' wavefunction in 
condensed matter, such as a bonding state or a lone- 
pair state, with a slight spatial extension. The suffix 
i of a Wannier state cj>i denotes its localization center. 
A bench mark is shown in the left panel of Fig. ^3 in 
which our calculations are compared with the conven- 
tional method for calculating eigenstates. The figure 
shows not only the result of the VR procedure but also 
that of another Wannier state method called 'perturba- 
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FIG. 9: Appearance of well-defined reconstructed surfaces 
after bending in cleavage path. The figures are drawn by 
projection from the [110] direction; (a): Bending from (111) 
plane to (110) plane. After the bending, the buckled (110) 
surface structures appears, as indicated by 'b'. (b)-(c): Two 
snapshots of the bending from (001) plane to (111) plane. The 
resultant (111) plane contains the 7r-bonded (2 x 1) structure 
with seven- and five- membered rings, indicated as '7' and '5' 
in (c). 



for decadesi&iSS and can be founded by the ab initio 
theory42i2ii2£ Consequently, group IV elements can be 
systematically described by a one-parameter energy scal- 
ing theoryii§iiL2LiS The scaling parameter, 'metallicity' 
a m , is defined as 



Here, (e p —e s ) is the energy difference between the p and s 
orbitals and t sp 3 is the transfer energy along a bulk (sp 3 ) 
bond (t sp s = \V ssa - 2^/W spa - W ppa \/A). A typical 
value for C is a m = 0.35^ The values for Si and Ge are 
similar, a m — 0.75 — 0.78 (see RefsiiSiii and references 
therein). 

Several surface structures among C, Si and Ge are sys- 
tematically explained by the energy scaling theory. An 
example is shown in the dimer geometry on (001) sur- 
faces. Ab initio calculations result in a symmetric dimer 
in the C case and similar asymmetric dimers in the Si 
and Ge cases (see Fig. 2 of Rcf. 44 , for example). In 
the right panel of Fig. ^| the above trend is reproduced 
with the present Hamiltonian by tuning the value of a m 
or (s p — e s ). Since the (lll)-(2 x 1) cleaved surface is 
commonly seen in Si and Gei^, the phenomena should 
hold a common mechanism in the present context of uni- 
versality. 

For validation of the present study, the surface en- 
ergy of several Si surfaces was calculated by the present 
Hamiltonian as follows and the ab initio values 1 ^ are 
given inside the parentheses; 7g ( ^ 1 2 = 1.58(1.41) J/m 2 , 
7 1 2 1 x 1 1 = 1.97 (1.44) J/m 2 and 7 j 3 1 u cklcd = 2.11 (1.70) J/m 2 , 
for the (001)-c(4x2) structure, the 7r-bonded (lll)-(2xl) 
structure and the buckled (110) structure, respectively. 
The present Hamiltonian reproduces the ab initio re- 
sults satisfactorily, particularly, the order of magnitude 
(7ooi 2 < 7m 1 < 7iio klcd )> though the absolute values are 
somewhat overestimated. The calculated surface energy 
of (001) surfaces also gives satisfactory results^ 



tive (PT) procedure '^LiSi The circle and square indicate 
the results of the conventional method and PT procedure, 
respectively, by a standard workstation^ The triangle 
and cross indicate the results of the VR procedure with 
32 CPUs and PT procedure with 512 CPUs, respectively, 
by a parallel computation system (SGI Origin 3800™). 
The parallelism is carried out by the OpenMP technique 
(www.openmp.org). All the results of our methods are 
'order- N\ or linearly proportional to the system size (N). 
The PT procedure is much faster than the VR procedure 
but its applicability is strictly limitedii§ii£ Several re- 
lated methods^ were also developed but are not used in 
this paper. 

Another methodological foundation is the transferable 
Hamiltonian 2 ^ in the Slater-Koster (tight-binding) form, 
applicable to various circumstances, e.g., crystals, liquid, 
defects and surfaces. Its success is based on the uni- 
versality of electronic structures, which has been known 



APPENDIX B: CONDITIONS OF CLEAVAGE 
SIMULATION 

This appendix explains the technical conditions for 
cleavage simulation. The time step of the molecular dy- 
namics was set at 3 fs. The center of gravity of the sample 
was fixed as a constraint. Sample surfaces were termi- 
nated by orientationally fixed sp 3 bonding states. This 
boundary condition corresponds to the situation in which 
the sample is embedded in a bulk (sp 3 -bonded) system. 
This situation is quite different from that with usual hy- 
drogen termination, in which the atomic structure is de- 
formed significantly from the tetrahedral geometry, ow- 
ing to the large deviation from the sp 3 bonding. 

Initial defect bonds were prepared for the cleavage 
'seed', which can be seen in Fig. |2Jb) or Fig. [HJb). Ac- 
tually, an additional short-range repulsive force was im- 
posed on the atom pairs of selected bond sites. The se- 
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FIG. 10: (Left) Computational time for molecular dynamics 
simulation with up to 11,315,021 Si atoms. Our calculations 
are compared with the conventional calculation for the eigen- 
states. See text for explanations. (Right) Dimer geometry on 
(001) surface of group IV elements (C, Si and Ge). The angle 
9 is the tilt angle and cj> is the angle between the surface dimer 
and the plane that contains the two back bonds of the upper 
atom. 



lected bond sites were, typically, two or three bond sites. 
Since the repulsive force is of short range, it will act on 
nothing after the bond breaking. 

Here, we discuss the way of imposing the external load 
on sample surfaces (boundaries). If the sample is suffi- 
ciently large, the cleavage phenomena will be determined 
by the intrinsic nature of the solid, not by the boundary 
condition of the sample. We have investigated various 
ways of imposing an external load^i and the discussion 
in this paper is based on the results that are not sensi- 
tive to the boundary conditions. We pick out the cases of 
Figs. |3 and IS] In Fig. [21 the external load is realized by 
the constrained movement. The atoms of the outermost 
two layers on the top and bottom sample surfaces were 
constrained under artificial movement in the [111] direc- 
tion. Since the cleavage propagation velocity is on the 
order of the Rayleigh wave velocity = 4.5nm/ps, the 
velocity of the artificial movement should be much slower 
than vr. If not, the atoms in the constraint motion (the 
atoms near the sample surfaces) are removed from the 
sample and no internal (cleavage) surface appears. In 
the case of Fig. [21 we chose vq = 0.16nm/ps, as a typ- 
ical value. In our simulations, the artificial movement 
was required to avoid forming multiple cleavage planes. 
The details are described elsewhereiSi Within the above 
requirements for the artificial movement, the conclusions 
of the present paper are not changed. The way of con- 
trolling the artificial movement is reflected on the shape 
of the sample boundary in Fig. [3 (a)-(e). 

Now we explain the case of Fig. [5] in which an external 
force field is imposed on selected atomic layers near the 
sample surfaces, except the left one. The thickness of 
the selected layers is 5 or 10 % of the sample length. The 



F(r) 



Kdl 



(Bl) 



with the relative coordinate r from the force center. 
The (two-dimensional) coordinate r is given on the (011) 
plane (See Fig. ^) . The force center (r = 0) is fixed at 
the position of the initial defect and is placed almost at 
the middle of the left sample surface. The length do is 
determined so that the volume per atom is given as dp 
(do ~ 3A). The force field is consistent with the concen- 
trated stress field given by the continuum mechanics of 
cleavage^ when its angular dependence is ignored. The 
angular dependence was ignored in our simulations, be- 
cause it is defined only for a planar crack and is not valid 
in the stepped or bent cleavage paths discussed in this pa- 
per. In Fig. [8] K = O-SMPay 7 ]!!, which is slightly larger 
than the critical value (if pr0 p = O-TMPa-^/m) estimated 
in Appendix In several other simulations, the field 
center was under artificial movement, owing to the fact 
that the crack tip moves with the cleavage propagation. 
We did not find, however, a systematic difference in the 
resultant cleavage behavior. 

Finally, in Figs. [21 and [SI the total kinetic energy was 
controlled by the thermostat method 46 with a typical 
temperature parameter of T=800K. The temperature pa- 
rameter T, however, does not correspond to experimen- 
tal temperature, since the simulation is a nonequilibrium 
process. In Fig. [21 for example, the constraint movement 
is introduced to the sample surfaces with the character- 
istic velocity vq- The temperature parameter T should 
be sufficiently large to propagate the introduced defor- 
mation into the internal region. 



APPENDIX C: ESTIMATION OF CRITICAL 
STRESS INTENSITY FACTOR 

The critical stress intensity factor for cleavage propa- 
gation, if prop, was estimated for the (lll)-(2x 1) cleavage 
mode with a smaller sample of 416 atoms. In the con- 
tinuum theory with linear elasticity^ the divergent stress 
field is presented at the crack tip. The center of the di- 
vergent field R c is given and the stress intensity factor 
K is introduced as the amplitude of the divergent field. 
Its critical value K prop is a measurement of the critical 
external load and was determined in several experimen- 
tal techniques and in the electronic structure calculations 
with hundreds of atoms. 

The present estimation of the critical factor K prop 
was carried out within the following quasi-static picture; 
when the field center R c is shifted by one atomic incre- 
ment, the crack tip should be shifted, also by one atomic 
increment, with the elementary bond-breaking and re- 
construction process. If the factor K is not sufficiently 
large (K < K prop ), however, the crack tip is not shifted. 
In practical simulations, the field center R c was chosen 
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in the middle of a bond layer, which is shown as the non- 
stepped path of Fig. Old). As constraints, outer atoms 
were fixed with the linear elastic displacement of plane 
strain^, determined by the factor K and the position 
R c . Internal (movable) atoms are relaxed in the finite 
temperature of 300 K for 2.4 ps. After the relaxation 
of the internal atoms, the position R c is shifted by 1/10 
of one atomic increment and the relaxation is restarted. 
The resultant cleaved surface shows the buckled (2 x 1) 
structure, as in Fig. 2 of Ref. 9 . In our simulations, the 
increase of the reconstructed surface area by one atomic 
increment occurs within the typical time scale of 0.1 ps. 

The present calculations resulted in a typical value of 
K prop = O.TMPay^, which agrees satisfactorily with the 
previous theoretical and experimental values of 0.41 — 
1.24MPaVm listed in Table I of Ref. 9 Since the calcu- 
lated surface energy in Appendix^ is somewhat overes- 
timated, for example, by approximately 40 % in 7]^ , we 



speculate that the present value of K wop is an overesti- 
mated one, owing to the limited freedom of the present 
Hamiltonian4I 
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